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We have investigated the possible cause-and-effect relationship due to stress transfer between two 
earthquakes that occurred near Christchurch, New Zealand, in September 2010 and in February 201 1. The 
Mw 7.1 Darfield (Canterbury) event took place along a previously unrecognized fault. The Mw 6.3 
Christchurch earthquake, generated by a thrust fault, occurred approximately five months later, 6 km 
south-east of Christchurch's city center. We have first measured the surface displacement field to retrieve 
the geometries of the two seismic sources and the slip distribution. In order to assess whether the first 
earthquake increased the likelihood of occurrence of a second earthquake, we compute the Coulomb Failure 
Function (CFF). We find that the maximum CFF increase over the second fault plane is reached exactly 
around the hypocenter of the second earthquake. In this respect, we may conclude that the Darfield 
earthquake contributed to promote the rupture of the Christchurch fault. 

Earthquakes interact. The research in the last decades demonstrated that over major active faults or fault 
systems, where seismologists registered the occurrence of an earthquake, the probability of occurrence of a 
second shock increases or decreases according to stress changes 1 4 . Indeed, a mainshock perturbs the stress 
state in other sections of the same fault or in adjacent faults: this theory is known as Coulomb Stress Triggering 5 . 
The hypothesis is that once an earthquake occurs, the stress does not dissipate, but it propagates in the surround- 
ing area, where it may increase the probability of occurrence of further earthquakes. Several examples can be 
found in the literature: in some cases, two seismic events, Earthquake-1 and Earthquake-2, occur over a long time 
span 4 ; in other cases they are temporally very close 6 ". The towns of Darfield (Canterbury) and Christchurch have 
been hit by strong earthquakes within a time span of a few months 912 . Earthquake-1 took place on September 03, 
2010 (at 16:35:46 UTC, Mw 7.1), while Earthquake-2 occurred on February 21, 2011 (at 23:51:43 UTC, Mw 6.3), 
a few kilometers more to the East. Between the two events, more than 4.000 aftershocks were recorded over 
the Greendale fault and other secondary structures. Could the second event have been triggered by a stress 
redistribution induced by the first one? The answer is in terms of probability change. Indeed stress change can 
trigger a second event if fault conditions are close to failure, otherwise it will hasten the occurrence of future 
earthquakes. 

New Zealand is located across the margin of the Australian and Pacific plates in the southern pacific. Here the 
relative obliquely convergent plate motion varies from 30 mm/yr in the south of the country to about 50 mm/yr 
in the northern part 13 . Starting from the Middle-Late Cretaceous and the Late Cretaceous-Paleocene episodes of 
rifting, New Zealand separated from Australia and Antarctica 14 since the Late Paleocene, during which a new 
extensional plate boundary initiated about 45 Ma 15 . Around 25 Ma the processes of oblique compression across 
the plate boundary initiated 16 , with an increase of convergence rates since the late Miocene. 

The study area is located within the Canterbury-Chathamas platform, characterized by a 1000-km-long dextral 
transpressive zone, represented by the Alpine fault and the Malborough fault belt (Figure la) ". This fault 
accommodates the relative plate motion between the NW-dipping Hikurangi and SE-dipping Puysegur subdue - 
tion zones 18 . On the Alpine fault, representative strike-slip rates are 25-30 mm/yr 19, 20 and only geological data 
provide evidence of large-to-great earthquake occurrence, with recurrence times of hundreds of years. The slip 
rates observed on the Malborough fault belt are greater than 1 mm/yr; here many moderate-to-large earthquakes 
occurred in historic times, including the Mw 8.1-8.2 1855 Wairarapa and the Mw 7.8 1931 Hawkes Bay 
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Figure 1 | (A) Tectonic setting of southern New Zealand island 14 . Black lines indicate the main tectonic structures discussed in the text. Red rectangle 
represent the study area. (B) Aftershocks distribution between the two main earthquakes on the September 3, 2010 and on February 21, 201 1. Red lines 
represent the activated faults. (C) Unwrapped interferograms in LOS geometry. Red stars are the hypocenter location (from CMT) of the two events; Focal 
mechanisms are also shown. 



earthquakes. This strike-slip province is characterized by several 
fault strands, hundreds of kilometers in length, associated with folds 
and reverse faults. Very few data are available about the distribution 
of segment boundaries on the main strands and to constrain the 
timing of paleoseismic events. 

The Mw 7.1 Darfield (Canterbury) earthquake occurred along 
a previously unrecognized east-west fault line, the strike-slip 
Greendale fault (Figure lb). In this area no active faults had pre- 
viously been mapped, nor are large historical earthquakes known to 
have occurred. The only active faults known in this region are located 
further west, at the foothills of the Southern Alps, where several 
M > 6-7 earthquakes have occurred in the past 150 years. The 



teleseismic moment tensor (USGS) and finite fault solutions 21,22 , 
providing a far-field observation of the earthquake, have indicated 
a dextral strike-slip fault (Figure lc) in agreement with both the 
orientation and the sense of slip of the documented surface rupture 11 . 
On the other hand, the near field seismological observations (first 
motion and regional moment tensor, 10 ) show a large reverse faulting 
component, in contrast with the teleseismic solutions. This is due to 
differences in the two measurement techniques, which analyze high 
(near field) and low (teleseismic) frequency waves and are therefore 
sensitive to small or large features respectively. A possible interpreta- 
tion is that the earthquake started as a thrust event in a smaller 
scale structure and continued over the main structure (Greendale 
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fault) with a strike-slip mechanism, as indicated by the teleseismic 
solutions, thus accommodating the regional stress and releasing the 
largest energy fraction. 

The 6.3 Mw Christchurch earthquake occurred on February 21 sl , 
2011, approximately 5 months after the Darfield earthquake. The 
hypocentre was approximately 6 km south-east of Christchurch's 
city center, at a depth of 5-6 km, generated by a blind ENE 
oblique-thrust, faulting at the easternmost limit of the Darfield after- 
shocks (Figure lb) 12 . No specific structure for this event is directly 
linked to the main fault of the 2010 mainshock. The focal mechanism 
solution indicated a right-lateral oblique thrust faulting mechanism 
(Figure lc). 

Roughly 20 years ago the Earth Sciences received the impact of the 
newborn Differential Synthetic Aperture Radar Interferometry 
(DInSAR) technique 23 . DInSAR has become a key element in 
multidisciplinary studies of earthquake 24 . In order to investigate 
Earthquake- 1 and Earthquake-2 we apply DInSAR, using two pairs 
of satellite images acquired by the Japanese mission ALOS 
(Advanced Land Observing Satellite), and its onboard PALSAR 
(Phased Array type L-band Synthetic Aperture Radar) sensor, along 
two adjacent tracks. We use the Shuttle Radar Topographic Mission 
digital elevation model 25 to remove topographic fringes from the 
interferograms (Figure lc). 

Results 

The map of the surface coseismic displacement (Supplementary 
Figure SI) of Earthquake- 1 shows a complex pattern, meaning that 
the Greendale fault and other secondary buried faults moved during 
the seismic event. Although this is not an evidence of the triggering of 
Earthquake-2, it can be noted that the surface displacement field 
extends near to the epicenter of the February 21 st earthquake. The 
latter occurred about 50 km East in a less complex scenario; we 
argue, based on DInSAR, that only a single fault is rupturing during 
the shock 12 (Supplementary Figure SI). Looking at both interfero- 
grams it is plain that PALSAR is an effective tool to capture the 
surface displacement field in case of moderate and strong earth- 
quakes 26, 27 . We apply an adaptive filter 28 in order to reduce the noise, 
and a Minimum Cost Flow (MCF) phase unwrapping algorithm 29 . 

Based on preliminary field observations 11 , the Greendale rupture 
zone is unequivocally the area of main deformation, with significant 
amounts of displacement. The average displacement is ~2.5m along 
the ~30km long main rupture and reaches a maximum of ~5m 
suggesting that the greatest energy fraction is released by the 
strike-slip Greendale fault. Smaller scale ruptures (blind thrusts) 
on nearby faults are also associated with this earthquake, but are 
expected to contribute less to the CFF estimation. 

The slip distribution for the Darfield earthquake is estimated by 
linearly inverting the sub sampled DInSAR deformation map. In a 
previous work' the authors have investigated the coseismic deforma- 
tion produced by the Darfield earthquake, following the main surface 
rupture, in order to constrain the Greendale fault and its related 
extensions. We fix our main fault in a similar way, following the 
documented surface rupture. We approximate the Greendale fault 
as a concatenation of three planar rectangular strike-slip fault planes, 
with 44km of total length and 12km of width. These are the main 
segment striking E-W and coincident with the major part of the 
mapped surface rupture, the NE-SW extension of the fault and 
the step-over (offset to the north) segment at the eastern end of 
the Greendale fault. The fault planes are subdivided into a discrete 
number of rectangular patches 30 . A Green's function matrix (kernel) 
is composed by imposing a unitary dislocation on each patch and 
subsequently collecting the E, N and vertical surface displacement 
components projected along the satellites line of sight (LOS). The 
linear inversion is performed adopting the Occam's smoothing 
scheme 31 , minimizing the chi-square and the second order derivative 
(Laplacian) to avoid large, unphysical oscillations in slip values. 



During the inversion process we also take into account three thrusts: 
the first one is near the hypocenter, the second at the NW end of the 
Greendale fault, distant from the main deformation field, and the 
third near the step-over segment. The largest part of the slip is con- 
centrated over the three segments of the main fault, with a peak value 
of 6.5m in the central one, suggesting that the majority of the energy 
is released by the Greendale fault (Figure 2). The resulting geodetic 
moment is ~5X 10 19 Nm (Mw —7.1) if a crustal rigidity of 30GPa is 
assumed in accordance with seismological estimates 32 . Our results 
for the slip distribution reveal a similar pattern to the one derived in 9 . 

In order to investigate the consistency of our results, we performed 
an uncertainty analysis, computing the covariance and resolution 
matrices for our model parameters 33 . The analysis was limited to 
the three main segments of the Greendale fault, for the sake of sim- 
plicity. The resolution matrix (R) can be used to estimate the spatial 
resolution at the location of any single fault patch in the model. If a 
single slip patch is perfectly resolved, the corresponding column in R 
will have a value of 1 at the main diagonal position, and the off 
diagonal elements will be zero. Results indicate (Supplementary 
Figure 3) that only the first row of patches is well resolved 
(R—0.9) and that the resolution decreases rapidly below 2km depth. 
The results for resolution indicate that only general slip features can 
be resolved at depth. The standard deviation is obtained from the 
square root of the elements on the main diagonal of the covariance 
matrix. The uncertainty bounds are acceptable and range between 0 
and 70cm everywhere in the fault plane, with the exception of the 
Western fault segment, where the error reaches higher values at 
depth (SE corner). 

The DInSAR technique offers a useful tool for fault characteriza- 
tion 34 . Indeed the fringe shape, rate and orientation can be related to 
fault parameters, such as geometric dimensions and orientation 
angles. We have used a novel approach 35 , based on the Okada 
model 30 and Neural Networks (NNs), to investigate the fault geo- 
metry of the Christchurch earthquake. One of the advantages of 
this method is that it rapidly achieves a determination of the 
rupture plane. 

Once the geometries of the two faults are defined (see Table 1), we 
focus our analysis to understand the role of the first earthquake in 
promoting the rupture of the second event through the evaluation of 
the Coulomb Failure Function (CFF). The CFF is obtained by com- 
puting the stress tensor corresponding to the elastic dislocation 
induced by the Canterbury earthquake, projecting it onto the rupture 
plane of the Christchurch earthquake, and evaluating the relative 
weights of normal and shear stresses, assuming an effective friction 
coefficient of 0.4. This value is in agreement with laboratory values of 
friction and moderate pore pressure when fluids are not fully 
expelled 4, 5| 36 . Positive or negative variations of the CFF indicate that 
the stress field is acting to promote, or oppose, the rupture, respec- 
tively. Our results (Fig. 3a) show that the rupture of the 2010 Darfield 
earthquake loaded a large portion of the crust with stress values 
exceeding 1 bar. If we take into account the three-dimensional loca- 
tion of the Christchurch rupture plane (Fig. 2), we see that that the 
shallower part of the fault (down to about 5km depth) has actually 
been unloaded by the Canterbury earthquake. On the contrary, the 
remaining, deeper part of the fault has been brought closer to rup- 
ture, with largest values of stress loading towards the south-western 
tip of the plane. The average CFF value on the loaded portion of the 
fault is over 0.01 MPa, with peak values exceeding 0.03 MPa. 
Remarkably, the peak stress increase occurs in the southwestern part 
of the fault, where rupture nucleation as occurred according to GNS 
localizations. These stress levels are definitely non-negligible, since a 
stress value of the order of 0.01 MPa is considered a threshold for 
effective triggering of seismic events 37 . We expect these estimates to 
be affected by several uncertainties. First, the geometric parameters 
of the reconstructed planes are known within a certain level of pre- 
cision, depending on the DinSAR data coverage and technical details 
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Figure 2 | (Top view) Selected frame of the East coast of South New Zealand. Surface projection of the fault planes adopted for the Darfield 2010 
(Earthquake- 1) and for the Christchurch 2011 (Earthquake-2) events. The sense of the slip for Earhquake-1 is right-lateral (180°) for the Greendale fault 
and its lateral extensions, as indicated by the red arrows, and thrust (90°) for the three minor faults. The black arrows indicate the dip direction of the three 
minor fault planes. Aftershocks distribution (black dots) corresponds to the time period between September 3 rd , 2010 and February 23" 1 , 2011 12 . (Bottom 
panel) A 3D perspective view for Earthquake- 1 and Earthquake 2. The largest part of the slip (max ~6.5m) is concentrated in the middle segment 
(Greendale fault) from 0 to 6km depth. Coulomb stress change is estimated for the Earthquake-2 fault plane. The red and black stars indicate the 
hypocenter of Earthquake- 1 and Earthquake-2 respectively (GNS Science). Both panels are in UTM WGS84 coordinate system. 
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Table 1 | The parameter retrieval of the Earthquake- 1 and Earthquake-2, obtained by means of Okada model and NNs. 
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of the modeling procedure. Also the assumption of an elastic half- 
space is likely to introduce a bias, given the tectonic framework of the 
region. Finally, postseismic deformations add a time-dependent 
stress perturbation that contributes to the total effect on the fault. 
Quantifying the impact of all these contributions turns out to be a 
quite complex task; however, they are not expected to alter the 
qualitative aspect of our conclusions. In Fig. 3b we compare the 
location of aftershocks with CFF variations on optimally oriented 
fault planes, assuming a background regional stress corresponding to 
an E-W compressional tectonics 38, 39 . Even if a complete analysis of 
the sequence seismicity would be beyond the scope of this work, we 
see from Fig. 3b that most of the aftershocks occur in loaded regions, 
suggesting the hypothesis that the sequence evolution may indeed be 
driven by stress redistribution mechanisms. 

Discussion 

The DInSAR results allowed us to investigate a relation between a 
pair of spatio-temporally close earthquakes. The interferograms pro- 
vided a first input to our analysis chain, in the form of surface dis- 
placement fields to be used for inversion modeling. Subsequently, 
using the fault geometries of both earthquakes and the slip distri- 
bution of the first one (Darfield earthquake, September 3, 2010), we 
calculated the CFF over the second fault plane. 

After the September 3 earthquake, a sequence began with after- 
shocks located along the Greendale fault and some hidden secondary 
faults. A magnitude 5.1 aftershock occurred on September 8 nearby 
the epicenter of the February 21 event. Even though this event 
occurred five months after the September 3 earthquake, some scien- 
tists consider the February 21 an aftershock caused by a fault rupture 
within the zone of aftershocks of the mainshock 40 . Research on earth- 
quake triggering investigated earthquake interaction at different 
scales: mainshock-mainshock and mainshock-aftershock ( 41 and 
references therein). Furthermore, it is now accepted that earthquake 
triggering occurs at all scales and that there is no mechanistic differ- 
ence between the origin of foreshocks, mainshocks and aftershocks. 
Based on such premises, Earthquake-2 can be interpreted as a main- 
shock along a second fault, promoted by the stress perturbation of 
Earthquake- 1, or as the largest aftershock of the sequence started 
with Earthquake- 1. It is however considered out of the scope of 
the present work to pursue this issue further, although it would be 



interesting from the statistical point of view. With our analysis we 
cannot state whether Earthquake-2 was (or was not) extremely 
unlikely to occur 40 without the preceding Earthquake- 1; however, 
the outcome of our work is that Earthquake- 1 has loaded the 
Earthquake-2 fault, bringing it closer to failure. 

Methods 

We selected a dataset of images from the Japanese ALOS satellite. ALOS has onboard 
a Synthetic Aperture Radar sensor, PALSAR, which is an active microwave sensor 
using L-band frequency to achieve cloud-free and day-and-night land observation. 
ALOS PALSAR allows a Fine Beam Single (FBS) polarization mode that achieves a 
spatial resolution of 5 m in azimuth and 7 to 44 m in ground range, according to the 
incidence angle, and a Fine Beam Dual (FBD) polarization mode, which, compared to 
the FBS, has the same azimuth resolution and a halved range resolution (14 to 88 m). 

In this study we use two pairs of ALOS PALSAR images along two adjacent ground 
tracks. The image pair relative to the September 3 rd , 2010, earthquake is in FBD mode 
with an incidence angle of about 38°. The pre-seismic image dates August 13, 2010 
and the post-seismic September 25, 2010. The image pair used to study the second 
earthquake, February 21, 2011, is in FBS mode with an incidence angle of about 38°. 
In this case pre- and post-seismic data are taken on January 1, 201 1 and February 25, 
2011, respectively. 

We apply the DInSAR technique to detect coseismic ground deformation. The 
interferograms have been computed with a square pixel of about 28 m achieved by 
applying a 3 by 8 multi-look factor in slant-range and azimuth respectively for FBD, 
and a 6 by 8 factor for FBS, in order to improve the signal-to-noise ratio. Furthermore, 
the 90-m Shuttle Radar Topography Mission (SRTM) digital elevation model 25 has 
been used to simulate and remove the topographic contribution of the interferometric 
phase. Both interferograms maintain a good coherence, which allows capturing most 
of the coseismic deformation, also for the second earthquake where a lot of damages 
occurred in the city, confirming the effectiveness of the PALSAR L-band sensor in the 
case of strong events 27 . In order to mitigate phase noise we apply an adaptive filter 28 , 
while to retrieve the Line Of Sight (LOS) displacements, we use a minimum cost flow 
phase unwrapping algorithm 29 . 

After mapping the coseismic ground deformation we perform a linear inversion in 
order to retrieve the slip distribution over the Greendale fault. We take into account 
the documented surface rupture by representing the Greendale fault as a concat- 
enation of three smaller faults. The geometric features of the rectangular planes are as 
follow: -1- central segment (Greendale fault), E-W orientation, Length 20km, Width 
12km, dip 87° -2- western segment (NW-SE extension), NW-SE orientation, Length 
12km, Width 12km, dip 87° -3- eastern segment (step over toward north), E-W 
orientation, Length 12km, Width 12km, dip 87°, -4- thrust fault near step over eastern 
segment, Length 8km, Width 8km, dip 65°, -5- thrust fault near hypocenter, Length 
8km, Width 8km, dip 75' J , -6- thrust fault near the NW-SE Greendale fault extension, 
Length 8km, Width 8km, dip 65°. The rake vector is fixed at 180 '(right lateral 
faulting) and 90°(thrust faulting) for -1-, -2-, -3- and -4-, -5-, -6- respectively 
(Supplementary Figure S4). Subsequently the fault model is subdivided into square 
patches of 2km per side. We compose the Green's function matrix (kernel) by 
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Figure 3 | (a) Coulomb stress change projected onto the Earthquake-2 rupture geometry estimated on a horizontal plane at 6 km depth. A black rectangle 
marks the surface projection of the 2011 Christchurch earthquake (Earthquake-2) . Yellow stars mark the location of the two events (GNS Science) .(b) Coulomb 
stress change induced by Earthquake- 1 on optimally oriented fault planes. Circles mark the epicentral location of aftershocks from the GNS catalogue. 
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imposing a unitary dislocation on each one of the patches, calculating the E, N and 
vertical surface displacement components, and projecting along the satellites LOS. 

In order to retrieve the fault geometry of the Christchurch earthquake we use an 
innovative approach based on Neural Networks 35 , that calculates the surface dis- 
placement field due to a dislocation on the fault plane at depth assuming an elastic half 
space, and using the analytical solutions provided by Okada 30 . This method has been 
chosen to investigate this second earthquake because it provides the major advantage 
of allowing a fast development of a preliminary model for the seismic source and thus 
a fast determination of the rupture plane for the Christchurch event. Moreover, it 
avoids the requirement for phase unwrapping. 

Neural networks (NNs) have already been recognized as being a powerful tool for 
inversion procedure in remote sensing applications. They are composed of an 
ensemble of nonlinear computational elements (called neurons) connected by the so- 
called synapses, each characterized by a synaptic weight. Compared to other methods, 
the use of NNs is often effective, because they can simultaneously address nonlinear 
dependences and complex physical behavior with reduced computational efforts and 
without the need of any a priori information. 

We then investigate the role of the first earthquake in promoting the rupture of the 
second event through the evaluation of the Coulomb Failure Function (CFF). First we 
compute the elastic strain tensor corresponding to the dislocation field induced by the 
Darfield (Canterbury) earthquake, using the analytical solutions provided by Okada 30 . 
By applying the relations of standard elasticity theory, the spatially- variable strain is 
converted into an incremental stress tensor that acts as a perturbation to the pre- 
existing (unknown) stress state of the crust. Once the perturbation to the stress tensor is 
known, its effect on a given fault mechanism can be assessed by evaluating the CFF 
variation, defined as ACFF — Ai + u(Arj n + AP), where At and Arj n are, respectively, 
the shear and normal stress changes, jj. is the friction coefficient and AP is the pore 
pressure change 37 . It is common to rewrite the definition of the CFF variation as ACFF 
= At +(i'Aa n , where [i is an effective friction coefficient that takes into account static 
friction, hydrostatic pressure and pore fluid pressure 4, 5 ; in our work we assumed 
u — 0.4. The knowledge of ACFF allows to establish whether the imposed stress field acts 
to promote (ACFF>0) or to oppose (ACFF<0) the dislocation on the considered fault. 
As a rule of thumb, in literature a CFF increase of 0.01 MPa (the magnitude of tidal 
loading) is considered a threshold for the effectiveness of triggering 4,s . 
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